Environmental indicator analysis-Tilefish
# Loading packages and data
library(ggplot2)
library(ecodata)
library(lubridate)
library(dplyr)
library(stringr)
library(marmap) # bathymetry
library(RColorBrewer)
library(ggnewscale)
library(sf)
library(cowplot)
library(tidyverse)
library(ggpubr)
library(sf)
library(ggdist)
library(ggpubr)
library(wesanderson)
library(raster)
library(mgcv)
#library(glmTMB)
library(ggpmisc)
library(mgcViz)
library(gratia)
# CPUE data (no env covariates)
gt_data_model_cpue <- read.csv(here::here('data/catch_data/gt_data_model_cpue.csv'))
#names(gt_data_model_cpue) <- tolower(names(gt_data_model_cpue))
# Add in column with cpue
# note: Paul indicated to use small mesh
gt_data_model_cpue <- gt_data_model_cpue %>%
rename_all(., .funs = tolower) %>%
mutate(mesh_bin = case_when(mesh_size <= 5.6 ~ 'SM',mesh_size >= 5.6 ~ 'LG',
TRUE ~ 'NA')) %>%
mutate(cpue_hr = sum_gt_catch/effort_dur)
# Catch data:
sfobs <-readRDS(here::here('data/catch_data/gold_tile_sf_ob_v1_temp_price.rds'))
sfob.env <- sfobs %>%
mutate(mesh_bin = case_when(mesh_size <= 5.6 ~ 'SM', mesh_size >= 5.6 ~ 'LG',
TRUE ~ 'NA'),
cpue_hr = SUM_GT_CATCH/effort_dur) %>%
filter(YEAR %in% c(1998:2022) & mesh_bin == 'SM') %>%
dplyr::select(DATE, YEAR, MONTH, YDAY,trip_id,hull_num, area, effort_dur,
SUM_GT_CATCH, cpue_hr, mesh_size, mesh_bin, depth, start_lat, start_lon,
bottomT, bottomT_avg, MIN_TEMP_C, MEAN_TEMP_C, MAX_TEMP_C,
TEMP_VARIANCE, TEMP_DEVIATION, MEAN_DPTH_M, tri, sed) %>%
mutate(YEAR = as.integer(YEAR)) %>%
rename_all(., .funs = tolower)
#write.csv(sfob.env, "C:/Users/stephanie.owen/Documents/tilefish_indicators/data/catch_data/sfob_env.csv", row.names=FALSE)
areas <- sort(unique(sfob.env$area))
catch.tally.ann <- sfob.env %>% # aggregate by year
group_by(year) %>%
summarise(ttl_sum = sum(sum_gt_catch))
# Length data from observer program
lengths <- read.csv(here::here('data/catch_data/gt_data_length_andy.csv'))
names(lengths) <- tolower(names(lengths))
# Recruitment estimates from 2021 report
recruit <- read.csv(here::here('data/assessment_data/tilefish_rec_estimate_2021.csv'))
# Merge SF/Obs catch data with recruit estimates:
catch_recruit <- cbind(recruit %>% filter(year %in% c(1998:2020)),
catch.tally.ann %>%
filter(year %in% c(1998:2020)) %>%
dplyr::select(ttl_sum))
# loading in shape files for maps
wd = here::here('shapefiles')
US.areas <- st_read(here::here('shapefiles/USA.shp'), quiet = TRUE)
canada.areas <- st_read(here::here('shapefiles/Canada.shp'), quiet = TRUE)
bts_strata <- st_read(here::here('shapefiles/NES_BOTTOM_TRAWL_STRATA.shp'),
quiet = TRUE)
# plot(bts_strata) # to see all bottom trawl strata
gtf_strata <- bts_strata %>%
filter(STRATUMA %in% c('01030', '01040', '01070', '01080', '01110', '01120',
'01140', '01150', '01670', '01680', '01710', '01720',
'01750', '01760')) # select just the gtf strata
# plot(gtf_strata)
bathy <- marmap::getNOAA.bathy(-81,-58, 27, 46)
bathy = fortify.bathy(bathy)Tilefish data
Catch data
Year-class strength is broadly defined as the number of fish spawned or hatched in a given year (Ricker, 1975).
Figure 1. Sum of catch (not accounting for effort), across years. Light blue shaded region represents the temporal range of observer records and red shaded region represents temporal range of study fleet records. The ‘purple’ region is where they overlap. Note that 2000-2005 for observer records had low sample size/number of vessels for tilefish, making the shaded region likely the best region to use for analysis. The vertical dashed lines represent strong year classes for this species (Nesslage et al. 2021). Red asterisk marks year that stock was deemed ‘re-built’.
# tot_catch == total (sum_catch) across hauls. so if tallying up annually,
# use sum_catch
# Strong year-classes: 1970, 1973, 1993, 1999, 2005, 2013
ggplot(catch.tally.ann, aes(x = factor(year), y = ttl_sum, group = 1))+
geom_rect(aes(xmin = '2007', xmax = '2022', ymin = -Inf, ymax = Inf),
fill = 'red', alpha = 0.02) +
geom_rect(aes(xmin = '2000', xmax = '2022', ymin = -Inf, ymax = Inf),
fill = 'lightblue', alpha = 0.05) +
geom_vline(xintercept = c('1993','1999', '2005', '2013'), lty = 2) +
geom_line(color = 'black', size = 1.5) +
annotate("text", label = "*",
x = 26, y = 14000, size = 8, colour = "red" )+
xlab('Year') +
ylab('Total sum tilefish catch') +
# facet_wrap(~month)+
theme(axis.text.x = element_text(color = 'black',
size = 12, angle = 45, vjust = 1, hjust=1)) +
ecodata::theme_facet()
CPUE
Figure 2. Catch-per-unit-effot for undirected trawl trips from the Study fleet and observer program. Zeros have been added using species association methodology (via jaccard index).
gt_data_model_cpue %>%
filter(mesh_bin == 'SM') %>% # note: Paul indicated to use small mesh
group_by(year, source) %>%
summarise(mean_cpue = mean(cpue_hr),.groups = 'drop') %>%
ggplot(aes(x=year,y=mean_cpue)) +
geom_line(lwd = 1) +
facet_wrap(~source) +
theme_bw()
gt_data_model_cpue %>%
filter(mesh_bin == 'SM') %>%
group_by(year) %>%
summarise(mean_cpue = mean(cpue_hr),.groups = 'drop') %>%
ggplot(aes(x=year,y=mean_cpue)) +
geom_line(lwd = 1) +
labs(title = 'Study fleet + Observer combined') +
theme_bw()
Maps (catch)
Tilefish catch locations (study fleet/observer)
yrs = sort(unique(gt_data_model_cpue$year))
#for(i in 1:length(yrs)){
yrmap <- function(yrs){
gt_data_model_cpue %>%
filter(start_lat < 42.5 & depth_est > 50 & year == yrs) %>%
mutate(bin = cut(year, seq(min(year), max(year) + 4, 4), right = FALSE)) %>%
ggplot() +
geom_sf(data = US.areas %>% st_as_sf(),color = 'gray20', fill = '#cbdbcd') +
geom_contour(data = bathy,
aes(x=x,y=y,z=-1*z),
breaks=c(50,100,150,200, Inf),
size=c(0.3),
col = 'darkgrey') +
stat_summary_2d(aes(x=start_lon, y=start_lat, z = cpue_hr),
binwidth=c(0.16666,0.16666)) +
scale_fill_viridis_c() +
theme(legend.position = "bottom",
legend.key.size = unit(0.2, "cm"),
legend.key.width = unit(1, "cm")) +
coord_sf(xlim = c(-75,-65.5), ylim = c(36,44), datum = sf::st_crs(4326)) +
labs(x = '', y = '', fill = 'CPUE') +
theme_bw()
}
for(i in 1:length(yrs)){
cat("\n#####", as.character(yrs[i]),"\n")
print(yrmap(yrs[i]))
cat("\n")
}2000
2001
2002
2003
2004
2005
2006
2007
2008
2009
2010
2011
2012
2013
2014
2015
2016
2017
2018
2019
2020
2021
2022
Recruitment estimate
Figure 3. Age-1 recruitment estimate from the 2021 tilefish assessment across all years
ggplot(recruit, aes(x = factor(year), y = recruit_est, group = 1))+
geom_rect(aes(xmin = '2007', xmax = '2022', ymin = -Inf, ymax = Inf),
fill = 'red', alpha = 0.02) +
geom_rect(aes(xmin = '2000', xmax = '2022', ymin = -Inf, ymax = Inf),
fill = 'lightblue', alpha = 0.05) +
geom_vline(xintercept = c('1993','1999', '2005', '2013'), lty = 2) +
geom_line(color = 'black', size = 1.5) +
annotate("text", label = "*",
x = 26, y = 14000, size = 8, colour = "red" )+
xlab('Year') +
ylab('Total sum tilefish catch') +
# facet_wrap(~month)+
theme(axis.text.x = element_text(color = 'black',
size = 12, angle = 45, vjust = 1, hjust=1)) +
ecodata::theme_facet()
### Thought: Should we isolate years associated w/strong year classes (or bad)
### for correlations and analyses? Figure 4. Recruitment estimates in focus years
- Paul has suggested that data used to calculate recruitment estimates was strongest from 2000 to present.
- Extended to 1998 to capture both the strong year class and El nino associated with that year.
ggplot(recruit %>% filter(year %in% c(1998:2022)),
aes(x = factor(year), y = recruit_est, group = 1))+
geom_vline(xintercept = c('1993','1999', '2005', '2013'), lty = 2) +
geom_line(color = 'black', size = 1.2) +
xlab('Year') +
ylab('Recruit estimates') +
theme(axis.text.x = element_text(color = 'black',
size = 12, angle = 45, vjust = 1, hjust=1)) +
ecodata::theme_facet()
Figure 5. Recruitment estimates and Study Fleet and Observer catch data. Black line denotes recruitment estimate, yellow denotes sum of annual catch data across both Study fleet and Observer programs.
options(scipen=999)
ggplot(catch_recruit) +
geom_line(aes(x = factor(year), y = recruit_est, group = 1),
col = 'black', size = 1.2) +
geom_line(aes(x = factor(year), y = ttl_sum*1000), size = 1.2,
color = 'goldenrod1', group = 1) +
scale_y_continuous(sec.axis = sec_axis(~./1000, name = 'Catch (lbs)')) +
geom_vline(xintercept = c('1993','1999', '2005', '2013'), lty = 2) +
xlab('Year') +
ylab('Recruit estimates') +
theme(axis.text.x = element_text(color = 'black',
size = 12, angle = 45, vjust = 1, hjust=1)) +
ecodata::theme_facet()
Length data
All lengths
Figure 6. Distribution of lengths Figure 7. Length frequencies Figure 8. Frequency of smaller individuals
# Define category breaks
size_breaks <- c(0,10,20,30,40, 50, 60, 70, 80, 90, 100)
# Making a function to bin the catches
label_interval <- function(breaks) {
paste0("(", breaks[1:length(breaks) - 1], "-", breaks[2:length(breaks)], ")")
}
labels = label_interval(size_breaks)
# length freq. table
tab = table(cut(lengths$lenanml,
breaks = size_breaks,
labels = label_interval(size_breaks)))
## Plot full distribution
ggplot(lengths,
aes(x = lenanml)) +
geom_bar(position = position_dodge(),
alpha = 0.4, fill= 'blue', color="black") +
xlab('Tilefish length (cm)') +
theme_bw() +
theme_facet()
# Plot length frequencies
barplot(tab, xlab = 'Length bins (mm)', main = '')
# Just the little ones
barplot(tab[1:3], xlab = 'Length bins (mm)', main = '')
Juveniles
Young of year - year 1 and 2 size class
ggplot(lengths %>% filter(lenanml <= 26),
aes(x = lenanml)) +
geom_bar(position = position_dodge(),
fill= 'slateblue', color="black") +
xlab('Tilefish length (cm)') +
theme_bw() +
theme_facet()
ggplot(lengths %>% filter(lenanml <= 26),
aes(x = lenanml, fill = numlen)) +
geom_bar(position = position_dodge(),
alpha = 0.4, fill= 'blue', color="black") +
xlab('Tilefish length (cm)') +
theme_bw() +
facet_wrap(~year) +
theme_facet()
Environmental data
The strong year classes for Golden Tilefish were 1993, 1998, 2005, 2013. Some of the underlying oceanographic processes that may be related to recruitment may influence habitat, retention/displacement and food availablity. These are explored below.
Habitat
Tilefish occupy a very narrow band of habitat conditions. Therefore, temperature and salinity may be of interest.
SST
# SST
sst<-read.csv(here::here('data/sst/sst_ts_gtf_strata.csv'))Maps (SST)
SST analysis
# SST by month with lag
sst.lag <- sst %>%
mutate(mean_sst_lag2 = lag(weighted_mean_sst,2),
mean_sst_lag3 = lag(weighted_mean_sst,3),
mean_sst_lag6 = lag(weighted_mean_sst,6))
# filter sfob.env
sfob <- sfob.env %>% filter(depth > 50, cpue_hr > 0) %>%
group_by(year,month) %>%
summarise(mean_cpue = mean(cpue_hr))
# Join with sf/ob data
sst.sfob <- dplyr::full_join(sst.lag, sfob, by = join_by(year, month)) %>%
dplyr::select(year, month, mean_cpue, mean_sst, weighted_mean_sst, mean_sst_lag2, mean_sst_lag3, mean_sst_lag6) %>%
tidyr::drop_na()
## SST no lag
ggplot2::ggplot(sst.sfob,
aes(x=mean_cpue, y=weighted_mean_sst)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean SST') +
labs(title = 'Mean SST no lag') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## SST summer
ggplot2::ggplot(sst.sfob %>% filter(month %in% c(7,8,9)),
aes(x=mean_cpue, y=weighted_mean_sst)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean SST') +
labs(title = 'Summer Mean SST') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## SST spring
ggplot2::ggplot(sst.sfob %>% filter(month %in% c(4,5,6)),
aes(x=mean_cpue, y=weighted_mean_sst)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean SST') +
labs(title = 'Spring Mean SST') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## 2 month lag
ggplot2::ggplot(sst.sfob,
aes(x=mean_cpue, y=mean_sst_lag2)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean SST') +
labs(title = 'Mean SST 2 month lag') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## 3 month lag
ggplot2::ggplot(sst.sfob,
aes(x=mean_cpue, y=mean_sst_lag3)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean SST') +
labs(title = 'Mean SST 3 month lag') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## 6 month lag
ggplot2::ggplot(sst.sfob,
aes(x=mean_cpue, y=mean_sst_lag6)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean SST') +
labs(title = 'Mean SST 6 month lag') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
SST Correlation and GAM
# Correlation and plot
lm_sst<-lm(weighted_mean_sst ~ mean_cpue, data=sst.sfob)
summary(lm_sst)##
## Call:
## lm(formula = weighted_mean_sst ~ mean_cpue, data = sst.sfob)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.271 -4.142 -1.321 4.245 12.278
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.77170 0.61536 25.630 < 0.0000000000000002 ***
## mean_cpue -0.28762 0.07307 -3.936 0.000122 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.216 on 164 degrees of freedom
## Multiple R-squared: 0.08631, Adjusted R-squared: 0.08074
## F-statistic: 15.49 on 1 and 164 DF, p-value: 0.0001221plot(weighted_mean_sst ~ mean_cpue, data=sst.sfob, pch=1, col="dodgerblue") + abline(lm_sst)
## integer(0)cor.test(sst.sfob$mean_cpue, sst.sfob$weighted_mean_sst, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: sst.sfob$mean_cpue and sst.sfob$weighted_mean_sst
## t = -3.9361, df = 164, p-value = 0.0001221
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4270047 -0.1480980
## sample estimates:
## cor
## -0.293792ggscatter(sst.sfob, x = "mean_cpue", y = "weighted_mean_sst",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.coeff.args = list(method="pearson"),
add.params=list(color="dodgerblue", fill="lightgray"),
xlab = "Mean CPUE/hr", ylab = "Weighted Mean SST")
# GAM
sst_gam <- gam(mean_cpue ~ s(weighted_mean_sst), data=sst.sfob, method="REML")
#gam.check(sst_gam)
#summary(sst_gam)
draw(sst_gam, residuals=TRUE)
Bottom Temperature
Figure 1. GLORYS vs in-situ bottom temperatures from study fleet vessels.
Figure 2. Bottom temperature (C) across years. Blue dots are in-situ data, red dots are from GLORYS.
ggplot2::ggplot(sfob.env, aes(x=bottomt, y=mean_temp_c)) +
geom_point(color="blue", alpha=0.1)+
geom_abline(intercept = 0, slope = 1) +
xlab('Bottom Temp (SF)') +
ylab('Bottom Temp (GLORYS)') +
theme_bw() 
ggplot2::ggplot(sfob.env, aes(x=bottomt, y=year)) +
geom_point(color="blue", alpha=0.1) +
geom_point(data = sfob.env, aes(x=mean_temp_c, y=year),
color="red", alpha=0.1) +
xlab('Bottom Temp') +
ylab('Year') +
labs(color = 'Source') +
theme_bw() 
Maps (Bottom T)
jet.colors <-colorRampPalette(c("blue", "#007FFF", "cyan","#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))
# select just years with study fleet bottom temps
sf.bt <- sfob.env %>% filter(year>2006 & depth > 50)
yrs = sort(unique(sf.bt$year))
#for(i in 1:length(yrs)){
yrmap <- function(yrs){
sf.bt %>% filter(year == yrs) %>%
ggplot() +
geom_sf(data = US.areas %>% st_as_sf(),color = 'gray20', fill = '#cbdbcd') +
geom_contour(data = bathy,
aes(x=x,y=y,z=-1*z),
breaks=c(50,100,150,200, Inf),
size=c(0.3),
col = 'darkgrey') +
stat_summary_2d(aes(x=start_lon, y=start_lat, z = bottomt),
binwidth=c(0.16666,0.16666)) +
scale_fill_gradientn(colors = jet.colors(20)) +
coord_sf(xlim = c(-75,-65.5), ylim = c(36,44), datum = sf::st_crs(4326)) +
labs(x = '', y = '', fill = 'Bottom temperature (°C)') +
theme_bw()
}
for(i in 1:length(yrs)){
cat("\n######", as.character(yrs[i]),"\n")
print(yrmap(yrs[i]))
cat("\n")
}2007
2008
2009
2010
2011
2012
2013
2014
2015
2016
2017
2018
2019
2020
2021
2022
Bottom T analysis
The following figures compare in-situ bottom temperature from the study-fleet data set to the recruitment estimate.
- Note here temperatures are averaged across all depths > 50 for each month.
# Note here temperatures are averaged across all depths > 50 for each month.
# Create in-situ bottom temps by month w/lag
df.lag = sfob.env %>% filter(year > 2006 & depth > 50) %>%
group_by(year,month) %>%
summarise(mean_dpth = mean(depth),
mean_bt = mean(bottomt)) %>%
mutate(mean_bt_lag2 = lag(mean_bt,2),
mean_bt_lag3 = lag(mean_bt,3),
mean_bt_lag6 = lag(mean_bt, 6))
# Join in-situ bottom temps w/assessment recruitment estimate
df.join = dplyr::full_join(recruit, df.lag, by = join_by(year)) %>%
dplyr::select(year, month, recruit_est, mean_dpth,
mean_bt, mean_bt_lag2, mean_bt_lag3, mean_bt_lag6) %>%
tidyr::drop_na()
# See what months have data
sort(unique(df.join$month))## [1] 7 8 9 10 11 12hist(df.join$month) # will group into spring/summer fall/winter categories
## spring/summer bottom temp no lag
ggplot2::ggplot(df.join %>% filter(month %in% c(4,5,6,7,8)),
aes(x=recruit_est, y=mean_bt)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean bottom temp (°C)')+
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## spring/summer bottom temp no lag
ggplot2::ggplot(df.join %>% filter(month %in% c(9,10,11,12)),
aes(x=recruit_est, y=mean_bt)) +
geom_point(color= 'black')+
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean bottom temp (°C)')+
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## spring/summer bottom temp 2 month lag
ggplot2::ggplot(df.join %>% filter(month %in% c(4,5,6,7,8)),
aes(x = recruit_est, y = mean_bt_lag2)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean bottom temp (°C)')+
labs(title = 'Lag 2 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## spring/summer bottom temp 2 month lag
ggplot2::ggplot(df.join %>% filter(month %in% c(9,10,11,12)),
aes(x = recruit_est, y = mean_bt_lag2)) +
geom_point(color= 'black')+
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean bottom temp (°C)')+
labs(title = 'Lag 2 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## spring/summer bottom temp 3 month lag
ggplot2::ggplot(df.join %>% filter(month %in% c(4,5,6,7,8)),
aes(x = recruit_est, y = mean_bt_lag3)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean bottom temp (°C)')+
labs(title = 'Lag 3 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## spring/summer bottom temp 3 month lag
ggplot2::ggplot(df.join %>% filter(month %in% c(9,10,11,12)),
aes(x = recruit_est, y = mean_bt_lag3)) +
geom_point(color= 'black')+
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean bottom temp (°C)')+
labs(title = 'Lag 3 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## spring/summer bottom temp 6 month lag
ggplot2::ggplot(df.join %>% filter(month %in% c(4,5,6,7,8)),
aes(x = recruit_est, y = mean_bt_lag6)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean bottom temp (°C)')+
labs(title = 'Lag 6 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## spring/summer bottom temp 6 month lag
ggplot2::ggplot(df.join %>% filter(month %in% c(9,10,11,12)),
aes(x = recruit_est, y = mean_bt_lag6)) +
geom_point(color= 'black')+
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean bottom temp (°C)')+
labs(title = 'Lag 6 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
The following figures compare in-situ bottom temperature from the study-fleet data set to the Study fleet and Observer tilefish catch data.
- Note here temperatures are averaged across all depths > 50 for each month.
#filter by year, depth, cpue (includes min, max, and sd of bottomt by month)
df.lag.cpue = sfob.env %>% filter(year > 2006 & depth > 50 & cpue_hr > 0) %>%
group_by(year,month) %>%
summarise(mean_dpth = mean(depth),
mean_bt = mean(bottomt),
min_bt = min(bottomt),
max_bt = max(bottomt),
sd_bt = sd(bottomt),
mean_cpue = mean(cpue_hr)) %>%
mutate(mean_bt_lag2 = lag(mean_bt,2),
mean_bt_lag3 = lag(mean_bt,3),
mean_bt_lag6 = lag(mean_bt, 6))
#omit NA bottomt
df.lag.sfob <- df.lag.cpue[!(is.na(df.lag.cpue$mean_bt)), ]
# See what months have data
sort(unique(df.lag.sfob$month))## [1] 1 2 3 4 5 6 7 8 9 10 11 12hist(df.lag.sfob$month) #winter/spring and summer/fall; summer only season with any non-zero relationship, except fall with 6 mo lag
## Winter/spring bottom temp no lag
ggplot2::ggplot(df.lag.sfob %>% filter(month %in% c(1,2,3,4,5,6)),
aes(x=mean_cpue, y=mean_bt)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean bottom temp (°C)') +
labs(title = 'Winter/spring no lag') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Summer/fall bottom temp no lag
ggplot2::ggplot(df.lag.sfob %>% filter(month %in% c(7,8,9,10,11,12)),
aes(x=mean_cpue, y=mean_bt)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean bottom temp (°C)')+
labs(title = 'Summer/fall no lag') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Summer bottom temp no lag
ggplot2::ggplot(df.lag.sfob %>% filter(month %in% c(7,8,9)),
aes(x=mean_cpue, y=mean_bt)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean bottom temp (°C)')+
labs(title = 'Summer no lag') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Winter/spring bottom temp 2 month lag
ggplot2::ggplot(df.lag.sfob %>% filter(month %in% c(1,2,3,4,5,6)),
aes(x = mean_cpue, y = mean_bt_lag2)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean bottom temp (°C)')+
labs(title = 'Winter/spring lag 2 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Summer/fall bottom temp 2 month lag
ggplot2::ggplot(df.lag.sfob %>% filter(month %in% c(7,8,9,10,11,12)),
aes(x = mean_cpue, y = mean_bt_lag2)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean bottom temp (°C)')+
labs(title = 'Summer/fall lag 2 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Summer bottom temp 2 month lag
ggplot2::ggplot(df.lag.sfob %>% filter(month %in% c(7,8,9)),
aes(x = mean_cpue, y = mean_bt_lag2)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean bottom temp (°C)')+
labs(title = 'Summer lag 2 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Winter/spring bottom temp 3 month lag
ggplot2::ggplot(df.lag.sfob %>% filter(month %in% c(1,2,3,4,5,6)),
aes(x = mean_cpue, y = mean_bt_lag3)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean bottom temp (°C)')+
labs(title = 'Winter/spring lag 3 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Summer/fall bottom temp 3 month lag
ggplot2::ggplot(df.lag.sfob %>% filter(month %in% c(7,8,9,10,11,12)),
aes(x = mean_cpue, y = mean_bt_lag3)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean bottom temp (°C)')+
labs(title = 'Summer/fall lag 3 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Summer bottom temp 3 month lag
ggplot2::ggplot(df.lag.sfob %>% filter(month %in% c(7,8,9)),
aes(x = mean_cpue, y = mean_bt_lag3)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean bottom temp (°C)')+
labs(title = 'Summer lag 3 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Summer/fall bottom temp 6 month lag
# No 6 month lag data for winter/spring
ggplot2::ggplot(df.lag.sfob %>% filter(month %in% c(7,8,9,10,11,12)),
aes(x = mean_cpue, y = mean_bt_lag6)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean bottom temp (°C)')+
labs(title = 'Summer/fall lag 6 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Summer bottom temp 6 month lag
ggplot2::ggplot(df.lag.sfob %>% filter(month %in% c(7,8,9)),
aes(x = mean_cpue, y = mean_bt_lag6)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean bottom temp (°C)')+
labs(title = 'Summer lag 6 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Fall bottom temp 6 month lag
ggplot2::ggplot(df.lag.sfob %>% filter(month %in% c(10,11,12)),
aes(x = mean_cpue, y = mean_bt_lag6)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean bottom temp (°C)')+
labs(title = 'Fall lag 6 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
Corrleation and figures for sf/ob data and in-situ bottom temps
lm_bt<-lm(mean_bt ~ mean_cpue, data=df.lag.sfob)
summary(lm_bt)##
## Call:
## lm(formula = mean_bt ~ mean_cpue, data = df.lag.sfob)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.6286 -1.0257 0.1391 1.1202 3.2579
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.55520 0.24264 47.623 <0.0000000000000002 ***
## mean_cpue 0.03508 0.03154 1.112 0.268
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.72 on 122 degrees of freedom
## Multiple R-squared: 0.01004, Adjusted R-squared: 0.001924
## F-statistic: 1.237 on 1 and 122 DF, p-value: 0.2682#plot(mean_bt ~ mean_cpue, data=df.lag.sfob, pch=1, col="dodgerblue") + abline(lm_bt)
cor.test(df.lag.sfob$mean_cpue, df.lag.sfob$mean_bt, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: df.lag.sfob$mean_cpue and df.lag.sfob$mean_bt
## t = 1.1123, df = 122, p-value = 0.2682
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.07749342 0.27170887
## sample estimates:
## cor
## 0.1001921ggscatter(df.lag.sfob, x = "mean_cpue", y = "mean_bt",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.coeff.args = list(method="pearson",
label.x=20, label.y=7.5),
add.params=list(color="dodgerblue", fill="lightgray"),
xlab = "Mean CPUE/hr", ylab = "Mean bottom temp (°C)")
#ggqqplot(df.lag.sfob$mean_bt, ylab="Mean BT")
#ggqqplot(df.lag.sfob$mean_cpue, ylab="Mean CPUE")SFOB insitu bottom temp GAM
bt_gam <- gam(mean_cpue ~ s(mean_bt), data=df.lag.sfob, method="REML")
#gam.check(bt_gam)
#summary(bt_gam)
draw(bt_gam, residuals=TRUE)
Salinity
Here we explore salinity from the GLORYS reanalysis model at three different depths
- 55 meters (relevant to larvae, recruits)
- 110 meters (relevant recruits, juveniles)
- 220 meters (relevant to juveniles, adults)
# Salinity Maps (Salinity)
jet.colors <-colorRampPalette(c("blue", "#007FFF", "cyan","#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))
# b <- brick(here::here('data/salinity/dd_sal_55_2000_2009.tif'))
# b.00 <- b[,,,1:365]
# select just years with study fleet bottom temps
sf.bt <- sfob.env %>% filter(year>2006 & depth > 50)
yrs = sort(unique(sf.bt$year))
#for(i in 1:length(yrs)){
yrmap <- function(yrs){
sf.bt %>% filter(year == yrs) %>%
ggplot() +
geom_sf(data = US.areas %>% st_as_sf(),color = 'gray20', fill = '#cbdbcd') +
geom_contour(data = bathy,
aes(x=x,y=y,z=-1*z),
breaks=c(50,100,150,200, Inf),
size=c(0.3),
col = 'darkgrey') +
# stat_summary_2d(aes(x=start_lon, y=start_lat, z = bottomt),
# binwidth=c(0.16666,0.16666)) +
# scale_fill_gradientn(colors = jet.colors(20)) +
coord_sf(xlim = c(-75,-65.5), ylim = c(36,44), datum = sf::st_crs(4326)) +
labs(x = '', y = '', fill = 'Salinty') +
theme_bw()
}
for(i in 1:length(yrs)){
cat("\n######", as.character(yrs[i]),"\n")
print(yrmap(yrs[i]))
cat("\n")
}2007
2008
2009
2010
2011
2012
2013
2014
2015
2016
2017
2018
2019
2020
2021
2022
Salinity analysis
Currents
Cross-shelf processes may influence the retention or displacement of tilefish during early life history stages. These are explored below.
Shelf water volume
Shelf water volume: A measure of the volume of water bounded inshore of the shelf-slope front. In this analysis, shelf water is defined as all water having salinity <34 psu. The position of the shelf-slope front varies inter-annually with the higher shelf water values indicating the front being pushed further towards the shelf break.
high shv: front pushed towards sbf low shv: front pushed inshore (more slope water on shelf)
Hypothesis: Higher recruitment success correlated with years of higher shelf water volume in spring/summer. These months months may be particularly important as that is when spawning is occurring and the position of the sbf may influence the position of larvae (away from spawning grounds).
Additional variables in this dataset are shelf water temperature and salinity which may also be indicative of habitat conditions.
# Shelf water volume
shlfvol <- read.csv(here::here('data/shelf_water_volume/ShelfWaterVolume_BSB_update.csv'))
# wrangling date info, converting doy to date and month
yrs <- as.vector(nrow(shlfvol))
shlfvol$Year <- as.character(shlfvol$Year)
for (i in 1:nrow(shlfvol)){
yrs[i] <- strsplit(shlfvol$Year, ".", fixed = TRUE)[[i]][1]
}
shlfvol$year <- yrs
shlfvol <- shlfvol %>% mutate(date_= as.Date(Year.Day-1,
origin=paste0(year, "-01-01")),
month= strftime(date_, "%m"),
day=strftime(date_,"%d"),
year = as.integer(year),
month = as.numeric(month))
# Create shw vol by month w/lag
df.lag = shlfvol %>%
group_by(year,month) %>%
summarise(mean_t = mean(ShW.T),
mean_s = mean(ShW.S),
mean_v = mean(ShW.Vol)) %>%
mutate(mean_t_lag2 = lag(mean_t,2),
mean_t_lag3 = lag(mean_t,3),
mean_t_lag6 = lag(mean_t,6),
mean_s_lag2 = lag(mean_s,2),
mean_s_lag3 = lag(mean_s,3),
mean_s_lag6 = lag(mean_s,6),
mean_v_lag2 = lag(mean_v,2),
mean_v_lag3 = lag(mean_v,3),
mean_v_lag6 = lag(mean_v,6))
# Join in-situ bottom temps w/assessment recruitment estimate
df.join = dplyr::full_join(recruit, df.lag, by = join_by(year)) %>%
dplyr::select(year, month, recruit_est, mean_t, mean_s, mean_v,
mean_t_lag2, mean_s_lag2, mean_v_lag2,
mean_t_lag3, mean_s_lag3, mean_v_lag3) %>%
tidyr::drop_na()
# See what months have data
sort(unique(df.join$month))## [1] 7 8 9 10 11hist(df.join$month) # will group into spring/summer fall/winter categories
## Shelf water volume no lag
ggplot2::ggplot(df.join,
aes(x=recruit_est, y=mean_v)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean shelf water volume') +
labs(title = 'Shelf water volume') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Shelf water temperature no lag
ggplot2::ggplot(df.join,
aes(x=recruit_est, y=mean_t)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean shelf water temperature') +
labs(title = 'Shelf water temperature') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Shelf water salinity no lag
ggplot2::ggplot(df.join,
aes(x=recruit_est, y=mean_s)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean shelf water salinity') +
labs(title = 'Shelf water salinity') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
With lags
2 Months
## Shelf water volume 2 month lag
ggplot2::ggplot(df.join,
aes(x=recruit_est, y=mean_v_lag2)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean shelf water volume') +
labs(title = 'Shelf water volume - lag 2 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Shelf water temperature 2 month lag
ggplot2::ggplot(df.join,
aes(x=recruit_est, y=mean_t_lag2)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean shelf water temperature') +
labs(title = 'Shelf water temperature - lag 2 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Shelf water salinity 2 month lag
ggplot2::ggplot(df.join,
aes(x=recruit_est, y=mean_s_lag2)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean shelf water salinity') +
labs(title = 'Shelf water salinity - lag 2 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
3 Months
## Shelf water volume 3 month lag
ggplot2::ggplot(df.join,
aes(x=recruit_est, y=mean_v_lag3)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean shelf water volume') +
labs(title = 'Shelf water volume - lag 3 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Shelf water temperature 3 month lag
ggplot2::ggplot(df.join,
aes(x=recruit_est, y=mean_t_lag3)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean shelf water temperature') +
labs(title = 'Shelf water temperature - lag 3 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Shelf water salinity 3 month lag
ggplot2::ggplot(df.join,
aes(x=recruit_est, y=mean_s_lag3)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean shelf water salinity') +
labs(title = 'Shelf water salinity - lag 3 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
Annual
# Create shw vol by year w/lag
df.lag = shlfvol %>%
group_by(year) %>%
summarise(mean_t = mean(ShW.T),
mean_s = mean(ShW.S),
mean_v = mean(ShW.Vol)) %>%
mutate(mean_t_lag2 = lag(mean_t,2),
mean_t_lag3 = lag(mean_t,3),
mean_t_lag6 = lag(mean_t,6),
mean_s_lag2 = lag(mean_s,2),
mean_s_lag3 = lag(mean_s,3),
mean_s_lag6 = lag(mean_s,6),
mean_v_lag2 = lag(mean_v,2),
mean_v_lag3 = lag(mean_v,3),
mean_v_lag6 = lag(mean_v,6))
# Join in-situ bottom temps w/assessment recruitment estimate
df.join = dplyr::full_join(recruit, df.lag, by = join_by(year)) %>%
dplyr::select(year,recruit_est, mean_t, mean_s, mean_v,
mean_t_lag2, mean_s_lag2, mean_v_lag2,
mean_t_lag3, mean_s_lag3, mean_v_lag3,
mean_t_lag6, mean_s_lag6, mean_v_lag6) %>%
tidyr::drop_na()
## Shelf water vol
ggplot2::ggplot(df.join,
aes(x=year, y=mean_v)) +
geom_point(color = 'black') +
geom_line(color = 'black') +
xlab('Year') +
ylab('Mean shelf water volume') +
labs(title = 'Shelf water volume') +
theme_bw()
ggplot2::ggplot() +
geom_line(data = df.join, aes(x=year, y=mean_s), color = 'red') +
geom_line(data = df.join,aes(x=year, y=mean_t*1), color = 'blue') +
ylim(30.0,34.0) +
scale_y_continuous(name = 'Sh.Water Salinity',
sec.axis = sec_axis(~./1, name = 'Sh.Water Temperature')) +
xlab('Year') +
labs(title = 'Shelf water salinity/temperature') +
theme_bw()
ggplot2::ggplot() +
geom_line(data = df.join, aes(x=year, y=mean_s), color = 'red') +
xlab('Year') +
ylab('Mean shelf water salinity') +
labs(title = 'Shelf water salinity') +
theme_bw()
## Shelf water vol no lag
ggplot2::ggplot(df.join,
aes(x=recruit_est, y=mean_v)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean shelf water volume') +
labs(title = 'Shelf water volume') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Shelf water temperature no lag
ggplot2::ggplot(df.join,
aes(x=recruit_est, y=mean_t)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean shelf water temperature') +
labs(title = 'Shelf water temperature') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Shelf water salinity no lag
ggplot2::ggplot(df.join,
aes(x=recruit_est, y=mean_s)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean shelf water salinity') +
labs(title = 'Shelf water salinity') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
With lags
## Shelf water vol 6 yr lag
ggplot2::ggplot(df.join,
aes(x=recruit_est, y=mean_v_lag6)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean shelf water volume') +
labs(title = 'Shelf water volume - lag 6 years') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Shelf water temperature 6 yr lag
ggplot2::ggplot(df.join,
aes(x=recruit_est, y=mean_t_lag6)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean shelf water temperature') +
labs(title = 'Shelf water temperature - lag 6 years') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Shelf water salinity 6 yr lag
ggplot2::ggplot(df.join,
aes(x=recruit_est, y=mean_s_lag6)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Recruitment estimate') +
ylab('Mean shelf water salinity') +
labs(title = 'Shelf water salinity - lag 6 years') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
Shelf water volume with Study fleet/Observer data
# Create shw vol by month w/lag
df.lag = shlfvol %>%
group_by(year,month) %>%
summarise(mean_t = mean(ShW.T),
mean_s = mean(ShW.S),
mean_v = mean(ShW.Vol)) %>%
mutate(mean_t_lag2 = lag(mean_t,2),
mean_t_lag3 = lag(mean_t,3),
mean_s_lag2 = lag(mean_s,2),
mean_s_lag3 = lag(mean_s,3),
mean_v_lag2 = lag(mean_v,2),
mean_v_lag3 = lag(mean_v,3))
#filter sfob.env
shlfvol.sfob = sfob.env %>% filter(depth > 50, cpue_hr > 0) %>%
group_by(year,month) %>%
summarise(mean_dpth = mean(depth),
mean_bt = mean(bottomt),
mean_cpue = mean(cpue_hr))
# Join with sf/ob data
df.shlfvol.sfob = dplyr::full_join(shlfvol.sfob, df.lag, by = join_by(year, month)) %>%
dplyr::select(year, month, mean_cpue, mean_t, mean_s, mean_v,
mean_t_lag2, mean_s_lag2, mean_v_lag2,
mean_t_lag3, mean_s_lag3, mean_v_lag3) %>%
tidyr::drop_na()
#See what months have data
sort(unique(df.shlfvol.sfob$month))## [1] 7 8 9 10 11hist(df.shlfvol.sfob$month) #summer and fall
No lags
## Shelf water volume no lag
ggplot2::ggplot(df.shlfvol.sfob,
aes(x=mean_cpue, y=mean_v)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean shelf water volume') +
labs(title = 'Shelf water volume no lag') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
##Shelf water volume summer
ggplot2::ggplot(df.shlfvol.sfob %>% filter(month %in% c(7,8,9)),
aes(x=mean_cpue, y=mean_v)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean shelf water volume') +
labs(title = 'Shelf water volume summer') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
##Shelf water volume fall
ggplot2::ggplot(df.shlfvol.sfob %>% filter(month %in% c(10,11,12)),
aes(x=mean_cpue, y=mean_v)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean shelf water volume') +
labs(title = 'Shelf water volume fall') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Shelf water temp no lag
ggplot2::ggplot(df.shlfvol.sfob,
aes(x=mean_cpue, y=mean_t)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean shelf water temperature') +
labs(title = 'Shelf water temperature no lag') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Shelf water salinity no lag
ggplot2::ggplot(df.shlfvol.sfob,
aes(x=mean_cpue, y=mean_s)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean shelf water salinity') +
labs(title = 'Shelf water salinity no lag') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
ggplot2::ggplot(df.shlfvol.sfob %>% filter(year >2007 & month %in% c(7,8,9,10,11,12)),
aes(x=mean_cpue, y=mean_v)) +
geom_point() +
geom_smooth(method = "lm", formula = y ~ x, size = 1, se = FALSE,
aes(colour = 'Linear')) +
geom_smooth(method = "lm", formula = y ~ x + I(x^2),
size = 1, se = FALSE, aes(colour = 'Quadratic')) +
geom_smooth(method = "loess", formula = y ~ x,
size = 1, se = FALSE, aes(colour = 'Loess')) +
geom_smooth(method = "gam", formula = y ~ s(x),
size = 1, se = FALSE, aes(colour = 'Gam')) +
geom_smooth(method = "gam", formula = y ~ s(x, k = 3),
size = 1, se = FALSE, aes(colour = 'Gam2')) +
labs(title = 'SFOB x M.S.W Summer/Fall (2007:2021)') +
scale_color_manual(name='Model',
breaks=c('Linear', 'Quadratic', 'Loess', 'Gam', 'Gam2'),
values=c('Linear'='black', 'Quadratic'='blue',
'Loess'='red', 'Gam' = 'green',
'Gam2' = 'purple')) +
theme_bw()
ggplot2::ggplot(df.shlfvol.sfob %>% filter(year >2007 & month %in% c(7,8,9,10,11,12)),
aes(x=mean_cpue, y=mean_s)) +
geom_point() +
geom_smooth(method = "lm", formula = y ~ x, size = 1, se = FALSE,
aes(colour = 'Linear')) +
geom_smooth(method = "lm", formula = y ~ x + I(x^2),
size = 1, se = FALSE, aes(colour = 'Quadratic')) +
geom_smooth(method = "loess", formula = y ~ x,
size = 1, se = FALSE, aes(colour = 'Loess')) +
geom_smooth(method = "gam", formula = y ~ s(x),
size = 1, se = FALSE, aes(colour = 'Gam')) +
geom_smooth(method = "gam", formula = y ~ s(x, k = 3),
size = 1, se = FALSE, aes(colour = 'Gam2')) +
labs(title = 'SFOB x M.S.W Summer/Fall (2007:2021)') +
scale_color_manual(name='Model',
breaks=c('Linear', 'Quadratic', 'Loess', 'Gam', 'Gam2'),
values=c('Linear'='black', 'Quadratic'='blue',
'Loess'='red', 'Gam' = 'green',
'Gam2' = 'purple')) +
theme_bw()
ggplot2::ggplot(df.shlfvol.sfob %>% filter(year >2007 & month %in% c(7,8,9,10,11,12)),
aes(x=mean_cpue, y=mean_t)) +
geom_point() +
geom_smooth(method = "lm", formula = y ~ x, size = 1, se = FALSE,
aes(colour = 'Linear')) +
geom_smooth(method = "lm", formula = y ~ x + I(x^2),
size = 1, se = FALSE, aes(colour = 'Quadratic')) +
geom_smooth(method = "loess", formula = y ~ x,
size = 1, se = FALSE, aes(colour = 'Loess')) +
geom_smooth(method = "gam", formula = y ~ s(x),
size = 1, se = FALSE, aes(colour = 'Gam')) +
geom_smooth(method = "gam", formula = y ~ s(x, k = 3),
size = 1, se = FALSE, aes(colour = 'Gam2')) +
labs(title = 'SFOB x M.S.W Summer/Fall (2007:2021)') +
scale_color_manual(name='Model',
breaks=c('Linear', 'Quadratic', 'Loess', 'Gam', 'Gam2'),
values=c('Linear'='black', 'Quadratic'='blue',
'Loess'='red', 'Gam' = 'green',
'Gam2' = 'purple')) +
theme_bw()
2 month lags
## Shelf water volume 2 month lag
ggplot2::ggplot(df.shlfvol.sfob,
aes(x=mean_cpue, y=mean_v_lag2)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean shelf water volume') +
labs(title = 'Shelf water volume - lag 2 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
##Shelf water volume summer - 2 month lag
ggplot2::ggplot(df.shlfvol.sfob %>% filter(month %in% c(7,8,9)),
aes(x=mean_cpue, y=mean_v_lag2)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean shelf water volume') +
labs(title = 'Shelf water volume summer - 2 month lag') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
##Shelf water volume fall - 2 month lag
ggplot2::ggplot(df.shlfvol.sfob %>% filter(month %in% c(10,11,12)),
aes(x=mean_cpue, y=mean_v_lag2)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean shelf water volume') +
labs(title = 'Shelf water volume fall - 2 month lag') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Shelf water temperature 2 month lag
ggplot2::ggplot(df.shlfvol.sfob,
aes(x=mean_cpue, y=mean_t_lag2)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean shelf water temperature') +
labs(title = 'Shelf water temperature - lag 2 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Shelf water salinity 2 month lag
ggplot2::ggplot(df.shlfvol.sfob,
aes(x=mean_cpue, y=mean_s_lag2)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean shelf water salinity') +
labs(title = 'Shelf water salinity - lag 2 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
3 month lags
## Shelf water volume 3 month lag
ggplot2::ggplot(df.shlfvol.sfob,
aes(x=mean_cpue, y=mean_v_lag3)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean shelf water volume') +
labs(title = 'Shelf water volume - lag 3 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
##Shelf water volume summer - 3 month lag
ggplot2::ggplot(df.shlfvol.sfob %>% filter(month %in% c(7,8,9)),
aes(x=mean_cpue, y=mean_v_lag3)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean shelf water volume') +
labs(title = 'Shelf water volume summer - 3 month lag') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
##Shelf water volume fall - 3 month lag
ggplot2::ggplot(df.shlfvol.sfob %>% filter(month %in% c(10,11,12)),
aes(x=mean_cpue, y=mean_v_lag3)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean shelf water volume') +
labs(title = 'Shelf water volume fall - 3 month lag') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Shelf water temperature 3 month lag
ggplot2::ggplot(df.shlfvol.sfob,
aes(x=mean_cpue, y=mean_t_lag3)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean shelf water temperature') +
labs(title = 'Shelf water temperature - lag 3 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Shelf water salinity 3 month lag
ggplot2::ggplot(df.shlfvol.sfob,
aes(x=mean_cpue, y=mean_s_lag3)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Mean shelf water salinity') +
labs(title = 'Shelf water salinity - lag 3 months') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
Correlations between sf/ob data and shelf water volume/temp/salinity
lm_shlf<-lm(mean_v ~ mean_cpue, data=df.shlfvol.sfob)
summary(lm_shlf)##
## Call:
## lm(formula = mean_v ~ mean_cpue, data = df.shlfvol.sfob)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1481.6 -362.9 195.8 268.4 1373.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3598.75 274.41 13.114 0.000000000256 ***
## mean_cpue 47.17 42.34 1.114 0.281
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 744.2 on 17 degrees of freedom
## Multiple R-squared: 0.06805, Adjusted R-squared: 0.01323
## F-statistic: 1.241 on 1 and 17 DF, p-value: 0.2807plot(mean_v ~ mean_cpue, data=df.shlfvol.sfob, pch=1, col="dodgerblue") + abline(lm_shlf)
## integer(0)cor.test(df.shlfvol.sfob$mean_cpue, df.shlfvol.sfob$mean_v, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: df.shlfvol.sfob$mean_cpue and df.shlfvol.sfob$mean_v
## t = 1.1141, df = 17, p-value = 0.2807
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.2193327 0.6393225
## sample estimates:
## cor
## 0.2608647ggscatter(df.shlfvol.sfob, x = "mean_cpue", y = "mean_v",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.coeff.args = list(method="pearson"),
add.params=list(color="dodgerblue", fill="lightgray"),
xlab = "Mean CPUE/hr", ylab = "Mean shelf water volume", title="Shelf Water Volume and CPUE Correlation")
# shelf water temperature
lm_shlf_t<-lm(mean_t ~ mean_cpue, data=df.shlfvol.sfob)
summary(lm_shlf_t)##
## Call:
## lm(formula = mean_t ~ mean_cpue, data = df.shlfvol.sfob)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.8024 -0.7921 -0.3091 0.5551 2.5960
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.61789 0.49332 31.659 <0.0000000000000002 ***
## mean_cpue -0.10849 0.07611 -1.425 0.172
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.338 on 17 degrees of freedom
## Multiple R-squared: 0.1068, Adjusted R-squared: 0.05421
## F-statistic: 2.032 on 1 and 17 DF, p-value: 0.1721plot(mean_t ~ mean_cpue, data=df.shlfvol.sfob, pch=1, col="dodgerblue") + abline(lm_shlf_t)
## integer(0)cor.test(df.shlfvol.sfob$mean_cpue, df.shlfvol.sfob$mean_t, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: df.shlfvol.sfob$mean_cpue and df.shlfvol.sfob$mean_t
## t = -1.4254, df = 17, p-value = 0.1721
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6800244 0.1496892
## sample estimates:
## cor
## -0.3267349ggscatter(df.shlfvol.sfob, x = "mean_cpue", y = "mean_t",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.coeff.args = list(method="pearson"),
add.params=list(color="dodgerblue", fill="lightgray"),
xlab = "Mean CPUE/hr", ylab = "Mean shelf water temperature", title="Shelf Water Temperature and CPUE Correlation")
# shelf water salinity
lm_shlf_s<-lm(mean_s ~ mean_cpue, data=df.shlfvol.sfob)
summary(lm_shlf_s)##
## Call:
## lm(formula = mean_s ~ mean_cpue, data = df.shlfvol.sfob)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.42230 -0.08417 0.00136 0.09558 0.40936
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.93715 0.07099 463.94 <0.0000000000000002 ***
## mean_cpue -0.03045 0.01095 -2.78 0.0128 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1925 on 17 degrees of freedom
## Multiple R-squared: 0.3125, Adjusted R-squared: 0.2721
## F-statistic: 7.728 on 1 and 17 DF, p-value: 0.01283plot(mean_s ~ mean_cpue, data=df.shlfvol.sfob, pch=1, col="dodgerblue") + abline(lm_shlf_s)
## integer(0)cor.test(df.shlfvol.sfob$mean_cpue, df.shlfvol.sfob$mean_s, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: df.shlfvol.sfob$mean_cpue and df.shlfvol.sfob$mean_s
## t = -2.78, df = 17, p-value = 0.01283
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.8080659 -0.1405127
## sample estimates:
## cor
## -0.5590427ggscatter(df.shlfvol.sfob, x = "mean_cpue", y = "mean_s",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.coeff.args = list(method="pearson"),
add.params=list(color="dodgerblue", fill="lightgray"),
xlab = "Mean CPUE/hr", ylab = "Mean shelf water salinity", title="Shelf Water Salinity and CPUE Correlation")
Shelf water volume GAM
#shelf vol GAM
shlfvol_gam <- gam(mean_cpue ~ s(mean_v), data=df.shlfvol.sfob, method="REML")
#gam.check(shlfvol_gam)
#summary(shlfvol_gam)
draw(shlfvol_gam, residuals=TRUE)
Gulf-stream index
Gulf stream index was calculated based on method presented by Pérez-Hernández and Joyce (2014). The gulf stream index (GSI) is a measure of the degrees latitude above the average Gulf Stream position based on ocean temperature at 200m (15 C) depth between 55W to 75W.
Positive values indicate a the mean position of the GS is more Northernly, whereas negative values indicate a more Southernly position.
Test code for GSI/recruit year lags
gsi.m <- read.csv(here::here('data/gulf_stream_index/mm_gsi_1954_2022_chen.csv'))
df <- dplyr::full_join(recruit, gsi.m %>%
group_by(year) %>%
filter(month %in% c(3:8)) %>%
summarise(m.gsi = mean(GSI),
sd.gsi = sd(GSI),
max.gsi = max(GSI),
min.gsi = min(GSI)),
by = join_by(year)) %>%
mutate(mean_gsi_lag1 = lag(m.gsi,1),
mean_gsi_lag2 = lag(m.gsi,2),
mean_gsi_lag3 = lag(m.gsi,3),
mean_gsi_lag6 = lag(m.gsi,6)) %>%
mutate(pos = ifelse(m.gsi > 0, 'Northerly', 'Southerly'),
n.pos = ifelse(m.gsi > 0, 1, 0))
gsi_df <- df[-c(51:69),]
## No lag
lm_rec_gsi<-lm(m.gsi ~ recruit_est, data=gsi_df)
summary(lm_rec_gsi)##
## Call:
## lm(formula = m.gsi ~ recruit_est, data = gsi_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.37145 -0.55510 -0.02879 0.59257 1.29704
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.29742643206 0.23294946950 1.277 0.208
## recruit_est -0.00000002737 0.00000014146 -0.193 0.847
##
## Residual standard error: 0.7183 on 48 degrees of freedom
## Multiple R-squared: 0.0007791, Adjusted R-squared: -0.02004
## F-statistic: 0.03742 on 1 and 48 DF, p-value: 0.8474cor.test(gsi_df$recruit_est, gsi_df$m.gsi, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: gsi_df$recruit_est and gsi_df$m.gsi
## t = -0.19345, df = 48, p-value = 0.8474
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3038986 0.2523967
## sample estimates:
## cor
## -0.02791194ggscatter(gsi_df, x = "recruit_est", y = "m.gsi",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.coeff.args = list(method="pearson", label.x=10, label.y=-1.0),
add.params=list(color="dodgerblue", fill="lightgray"),
xlab = "Recruitment estimate", ylab = "Mean Yearly GSI")
#GSI GAM
gsi_gam <- gam(recruit_est ~ s(m.gsi), data=gsi_df, method="REML")
#gam.check(gsi_gam)
#summary(gsi_gam)
draw(gsi_gam, residuals=TRUE)
## 1 year lag
lm_rec_gsi1<-lm(mean_gsi_lag1 ~ recruit_est, data=gsi_df)
summary(lm_rec_gsi1)##
## Call:
## lm(formula = mean_gsi_lag1 ~ recruit_est, data = gsi_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.41297 -0.51227 -0.02784 0.51398 1.37455
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.5828359751 0.2265633341 2.573 0.0133 *
## recruit_est -0.0000002308 0.0000001370 -1.685 0.0986 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6948 on 47 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.05697, Adjusted R-squared: 0.03691
## F-statistic: 2.839 on 1 and 47 DF, p-value: 0.0986cor.test(gsi_df$recruit_est, gsi_df$mean_gsi_lag1, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: gsi_df$recruit_est and gsi_df$mean_gsi_lag1
## t = -1.6851, df = 47, p-value = 0.0986
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4871862 0.0455661
## sample estimates:
## cor
## -0.2386888ggscatter(gsi_df, x = "recruit_est", y = "mean_gsi_lag1",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.coeff.args = list(method="pearson", label.x=10, label.y=-1.0),
add.params=list(color="dodgerblue", fill="lightgray"),
xlab = "Recruitment estimate", ylab = "Mean Yearly GSI",
title = "GSI 1 year lag")
#GSI GAM
gsi_gam1 <- gam(recruit_est ~ s(mean_gsi_lag1), data=gsi_df, method="REML")
#gam.check(gsi_gam)
#summary(gsi_gam)
draw(gsi_gam1, residuals=TRUE)
## 2 year lag
lm_rec_gsi2<-lm(recruit_est ~ mean_gsi_lag2, data=gsi_df)
summary(lm_rec_gsi2)##
## Call:
## lm(formula = recruit_est ~ mean_gsi_lag2, data = gsi_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1085718 -545396 -139089 297537 2263014
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1550432 109190 14.199 <0.0000000000000002 ***
## mean_gsi_lag2 -271881 152993 -1.777 0.0822 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 722300 on 46 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.06424, Adjusted R-squared: 0.0439
## F-statistic: 3.158 on 1 and 46 DF, p-value: 0.08217cor.test(gsi_df$recruit_est, gsi_df$mean_gsi_lag2, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: gsi_df$recruit_est and gsi_df$mean_gsi_lag2
## t = -1.7771, df = 46, p-value = 0.08217
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.50148026 0.03305454
## sample estimates:
## cor
## -0.2534606ggscatter(gsi_df, x = "recruit_est", y = "mean_gsi_lag2",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.coeff.args = list(method="pearson", label.x=10, label.y=-1.0),
add.params=list(color="dodgerblue", fill="lightgray"),
xlab = "Recruitment estimate", ylab = "Mean Yearly GSI",
title = "GSI 2 year lag")
#GSI GAM
gsi_gam2 <- gam(recruit_est ~ s(mean_gsi_lag2), data=gsi_df, method="REML")
#gam.check(gsi_gam)
#summary(gsi_gam)
draw(gsi_gam2, residuals=TRUE)
## 3 year lag
lm_rec_gsi3<-lm(mean_gsi_lag3 ~ recruit_est, data=gsi_df)
summary(lm_rec_gsi3)##
## Call:
## lm(formula = mean_gsi_lag3 ~ recruit_est, data = gsi_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.26264 -0.49635 0.02971 0.46407 1.32520
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3431242573 0.2253764752 1.522 0.135
## recruit_est -0.0000001015 0.0000001351 -0.752 0.456
##
## Residual standard error: 0.6835 on 45 degrees of freedom
## (3 observations deleted due to missingness)
## Multiple R-squared: 0.01241, Adjusted R-squared: -0.009541
## F-statistic: 0.5653 on 1 and 45 DF, p-value: 0.4561cor.test(gsi_df$recruit_est, gsi_df$mean_gsi_lag3, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: gsi_df$recruit_est and gsi_df$mean_gsi_lag3
## t = -0.75185, df = 45, p-value = 0.4561
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3861955 0.1815940
## sample estimates:
## cor
## -0.1113811ggscatter(gsi_df, x = "recruit_est", y = "mean_gsi_lag3",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.coeff.args = list(method="pearson", label.x=10, label.y=-1.0),
add.params=list(color="dodgerblue", fill="lightgray"),
xlab = "Recruitment estimate", ylab = "Mean Yearly GSI",
title = "GSI 3 year lag")
#GSI GAM
gsi_gam3 <- gam(recruit_est ~ s(mean_gsi_lag3), data=gsi_df, method="REML")
#gam.check(gsi_gam)
#summary(gsi_gam)
draw(gsi_gam3, residuals=TRUE)
## 6 year lag
lm_rec_gsi6<-lm(mean_gsi_lag6 ~ recruit_est, data=gsi_df)
summary(lm_rec_gsi6)##
## Call:
## lm(formula = mean_gsi_lag6 ~ recruit_est, data = gsi_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.19710 -0.47456 0.00689 0.35780 1.30988
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2157713868 0.2100827803 -1.027 0.3103
## recruit_est 0.0000002303 0.0000001259 1.828 0.0746 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6331 on 42 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.07372, Adjusted R-squared: 0.05167
## F-statistic: 3.343 on 1 and 42 DF, p-value: 0.07461cor.test(gsi_df$recruit_est, gsi_df$mean_gsi_lag6, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: gsi_df$recruit_est and gsi_df$mean_gsi_lag6
## t = 1.8283, df = 42, p-value = 0.07461
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.0275832 0.5260005
## sample estimates:
## cor
## 0.2715206ggscatter(gsi_df, x = "recruit_est", y = "mean_gsi_lag6",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.coeff.args = list(method="pearson", label.x=10, label.y=-1.0),
add.params=list(color="dodgerblue", fill="lightgray"),
xlab = "Recruitment estimate", ylab = "Mean Yearly GSI",
title = "GSI 6 year lag")
#GSI GAM
gsi_gam6 <- gam(recruit_est ~ s(mean_gsi_lag6), data=gsi_df, method="REML")
#gam.check(gsi_gam)
#summary(gsi_gam)
draw(gsi_gam6, residuals=TRUE)
CPUE/Year lags
cpue_yr<-read.csv(here::here('data/catch_data/cpue_yr.csv'))
sfob_yr<- dplyr::full_join(cpue_yr, gsi.m %>%
group_by(year) %>%
filter(month %in% c(3:8)) %>%
summarise(m.gsi = mean(GSI),
sd.gsi = sd(GSI),
max.gsi = max(GSI),
min.gsi = min(GSI)),
by = join_by(year)) %>%
mutate(mean_gsi_lag1 = lag(m.gsi,1),
mean_gsi_lag2 = lag(m.gsi,2),
mean_gsi_lag3 = lag(m.gsi,3)) %>%
mutate(pos = ifelse(m.gsi > 0, 'Northerly', 'Southerly'),
n.pos = ifelse(m.gsi > 0, 1, 0))
gsi_yr <- sfob_yr[-c(26:69),]
## No lag
lm_cpue_gsi<-lm(m.gsi ~ mean_cpue_yr, data=gsi_yr)
summary(lm_cpue_gsi)##
## Call:
## lm(formula = m.gsi ~ mean_cpue_yr, data = gsi_yr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.4333 -0.4844 0.1177 0.5345 0.9829
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.55942 0.16989 3.293 0.00319 **
## mean_cpue_yr 0.02032 0.03260 0.623 0.53915
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7111 on 23 degrees of freedom
## Multiple R-squared: 0.01662, Adjusted R-squared: -0.02614
## F-statistic: 0.3886 on 1 and 23 DF, p-value: 0.5391cor.test(gsi_yr$mean_cpue_yr, gsi_yr$m.gsi, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: gsi_yr$mean_cpue_yr and gsi_yr$m.gsi
## t = 0.62341, df = 23, p-value = 0.5391
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.2805127 0.4986389
## sample estimates:
## cor
## 0.1289061ggscatter(gsi_yr, x = "mean_cpue_yr", y = "m.gsi",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.coeff.args = list(method="pearson", label.x=10, label.y=-1.0),
add.params=list(color="dodgerblue", fill="lightgray"),
xlab = "Mean Yearly CPUE", ylab = "Mean Yearly GSI",
title="GSI no lag")
#GSI GAM
gsi_cpue_gam <- gam(mean_cpue_yr ~ s(m.gsi), data=gsi_yr, method="REML")
#gam.check(gsi_gam)
#summary(gsi_gam)
draw(gsi_cpue_gam, residuals=TRUE)
## 1 year lag
lm_cpue_gsi1<-lm(mean_gsi_lag1 ~ mean_cpue_yr, data=gsi_yr)
summary(lm_cpue_gsi1)##
## Call:
## lm(formula = mean_gsi_lag1 ~ mean_cpue_yr, data = gsi_yr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.44875 -0.41955 0.04497 0.52003 0.95976
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.63929 0.16971 3.767 0.00106 **
## mean_cpue_yr -0.01995 0.03297 -0.605 0.55144
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7083 on 22 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.01636, Adjusted R-squared: -0.02835
## F-statistic: 0.3659 on 1 and 22 DF, p-value: 0.5514cor.test(gsi_yr$mean_cpue_yr, gsi_yr$mean_gsi_lag1, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: gsi_yr$mean_cpue_yr and gsi_yr$mean_gsi_lag1
## t = -0.60489, df = 22, p-value = 0.5514
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5052331 0.2904798
## sample estimates:
## cor
## -0.1279046ggscatter(gsi_yr, x = "mean_cpue_yr", y = "mean_gsi_lag1",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.coeff.args = list(method="pearson", label.x=10, label.y=-1.0),
add.params=list(color="dodgerblue", fill="lightgray"),
xlab = "Mean Yearly CPUE", ylab = "Mean Yearly GSI",
title="GSI 1 year lag")
#GSI GAM
gsi_cpue_gam1 <- gam(mean_cpue_yr ~ s(mean_gsi_lag1), data=gsi_yr, method="REML")
#gam.check(gsi_gam)
#summary(gsi_gam)
draw(gsi_cpue_gam1, residuals=TRUE)
## 2 year lag
lm_cpue_gsi2<-lm(mean_gsi_lag2 ~ mean_cpue_yr, data=gsi_yr)
summary(lm_cpue_gsi2)##
## Call:
## lm(formula = mean_gsi_lag2 ~ mean_cpue_yr, data = gsi_yr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.40538 -0.44659 -0.03594 0.56725 1.01038
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.540783 0.172542 3.134 0.00501 **
## mean_cpue_yr 0.004872 0.033200 0.147 0.88474
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7121 on 21 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.001024, Adjusted R-squared: -0.04655
## F-statistic: 0.02153 on 1 and 21 DF, p-value: 0.8847cor.test(gsi_yr$mean_cpue_yr, gsi_yr$mean_gsi_lag2, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: gsi_yr$mean_cpue_yr and gsi_yr$mean_gsi_lag2
## t = 0.14673, df = 21, p-value = 0.8847
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3852810 0.4384223
## sample estimates:
## cor
## 0.0320038ggscatter(gsi_yr, x = "mean_cpue_yr", y = "mean_gsi_lag2",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.coeff.args = list(method="pearson", label.x=10, label.y=-1.0),
add.params=list(color="dodgerblue", fill="lightgray"),
xlab = "Mean Yearly CPUE", ylab = "Mean Yearly GSI",
title="GSI 2 year lag")
#GSI GAM
gsi_cpue_gam2 <- gam(mean_cpue_yr ~ s(mean_gsi_lag2), data=gsi_yr, method="REML")
#gam.check(gsi_gam)
#summary(gsi_gam)
draw(gsi_cpue_gam2, residuals=TRUE)
## 3 year lag
lm_cpue_gsi3<-lm(mean_gsi_lag3 ~ mean_cpue_yr, data=gsi_yr)
summary(lm_cpue_gsi3)##
## Call:
## lm(formula = mean_gsi_lag3 ~ mean_cpue_yr, data = gsi_yr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.41731 -0.41127 -0.01401 0.49824 1.06394
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.46678 0.17512 2.665 0.0149 *
## mean_cpue_yr 0.02338 0.03312 0.706 0.4885
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7105 on 20 degrees of freedom
## (3 observations deleted due to missingness)
## Multiple R-squared: 0.0243, Adjusted R-squared: -0.02448
## F-statistic: 0.4982 on 1 and 20 DF, p-value: 0.4885cor.test(gsi_yr$mean_cpue_yr, gsi_yr$mean_gsi_lag3, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: gsi_yr$mean_cpue_yr and gsi_yr$mean_gsi_lag3
## t = 0.7058, df = 20, p-value = 0.4885
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.2844097 0.5418846
## sample estimates:
## cor
## 0.1558916ggscatter(gsi_yr, x = "mean_cpue_yr", y = "mean_gsi_lag3",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.coeff.args = list(method="pearson", label.x=10, label.y=-1.0),
add.params=list(color="dodgerblue", fill="lightgray"),
xlab = "Mean Yearly CPUE", ylab = "Mean Yearly GSI",
title="GSI 3 year lag")
#GSI GAM
gsi_cpue_gam3 <- gam(mean_cpue_yr ~ s(mean_gsi_lag3), data=gsi_yr, method="REML")
#gam.check(gsi_gam)
#summary(gsi_gam)
draw(gsi_cpue_gam3, residuals=TRUE)
#gsi.m <- read.csv(here::here('C:/Users/stephanie.owen/Documents/tilefish_indicators/gulf_stream_index/mm_gsi_1954_2022_chen.csv'))
#subset only GSI, month, year and 1997>
#gsi<-gsi.m[c("GSI","year","month")] %>%
# filter(year > 1997)
#write.csv(gsi, "C:/Users/stephanie.owen/Documents/tilefish_indicators/gulf_stream_index/gsi.csv", row.names=FALSE)
#filter sfob.env
#gsi.sfob <- sfob.env %>% filter(depth > 50) %>%
# group_by(year,month) %>%
#summarise(mean_cpue = mean(cpue_hr))
#write.csv(gsi.sfob, "C:/Users/stephanie.owen/Documents/tilefish_indicators/gulf_stream_index/gsi_sfob.csv", row.names=FALSE)
#new df wouldn't recognize month when joining the gsi and sfob data so I did it manually
gsi <- read.csv(here::here('C:/Users/stephanie.owen/Documents/tilefish_indicators/data/gulf_stream_index/gsi_sfob.csv'))
gsi_nozero <- filter(gsi, mean_cpue > 0)
gsi.sfob<-mutate(gsi_nozero, pos = ifelse(GSI > 0, 'Northerly', 'Southerly'),
n.pos = ifelse(GSI > 0, 1, 0))GSI by Month and Season
# GSI all months (1998-2022)
ggplot(gsi.sfob ,aes(x=mean_cpue, y=GSI)) +
geom_point(color = 'black') +
facet_wrap(~month)+
xlab('Mean CPUE/hr') +
ylab('Gulf stream position anomaly') +
labs(title = 'Gulf stream index by month') +
geom_hline(yintercept = 0, lty = 2) +
theme_bw()
#GSI Jan-Apr
filter(gsi.sfob, month %in% c(1:4)) %>%
ggplot(.,aes(x=mean_cpue, y=GSI)) +
geom_point(color = 'black') +
xlab('Mean CPUE/hr') +
ylab('Gulf stream position anomaly') +
labs(title = 'Gulf stream index (Jan:Mar)') +
geom_hline(yintercept = 0, lty = 2) +
theme_bw()
#GSI Sept-Dec
filter(gsi.sfob, month %in% c(9:12)) %>%
ggplot(.,aes(x=mean_cpue, y=GSI)) +
geom_point(color = 'black') +
xlab('Mean CPUE/hr') +
ylab('Gulf stream position anomaly') +
labs(title = 'Gulf stream index (Apr:Jul)') +
geom_hline(yintercept = 0, lty = 2) +
theme_bw()
# All months
ggplot(gsi.sfob, aes(x=mean_cpue, y=GSI)) +
geom_point(color = 'black') +
geom_hline(yintercept = 0, lty = 2)+
labs(title = 'All months') +
xlab('Mean CPUE/hr') +
ylab('Gulf stream position anomaly') +
geom_smooth(method = "lm", formula = y ~ x, size = 1, se = FALSE,
aes(colour = 'Linear')) +
geom_smooth(method = "lm", formula = y ~ x + I(x^2),
size = 1, se = FALSE, aes(colour = 'Quadratic')) +
geom_smooth(method = "loess", formula = y ~ x,
size = 1, se = FALSE, aes(colour = 'Loess')) +
geom_smooth(method = "gam", formula = y ~ s(x),
size = 1, se = FALSE, aes(colour = 'Gam')) +
geom_smooth(method = "gam", formula = y ~ s(x, k = 3),
size = 1, se = FALSE, aes(colour = 'Gam2')) +
scale_color_manual(name='Model',
breaks=c('Linear', 'Quadratic', 'Loess', 'Gam', 'Gam2'),
values=c('Linear'='black', 'Quadratic'='blue',
'Loess'='red', 'Gam' = 'green',
'Gam2' = 'purple')) +
theme_bw()
#ANOVA
pos_aov_cpue <- aov(mean_cpue ~ month * pos,
data = gsi.sfob)
# look at effects and interactions
summary(pos_aov_cpue)## Df Sum Sq Mean Sq F value Pr(>F)
## month 1 61.1 61.14 5.608 0.0192 *
## pos 1 27.0 27.00 2.476 0.1178
## month:pos 1 1.1 1.11 0.102 0.7498
## Residuals 143 1559.1 10.90
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tidy_pos_aov_cpue <- broom::tidy(pos_aov_cpue)
tidy_pos_aov_cpue## # A tibble: 4 × 6
## term df sumsq meansq statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 month 1 61.1 61.1 5.61 0.0192
## 2 pos 1 27.0 27.0 2.48 0.118
## 3 month:pos 1 1.11 1.11 0.102 0.750
## 4 Residuals 143 1559. 10.9 NA NAggplot(data = gsi.sfob,
aes(x = mean_cpue, y = GSI, fill = pos,
group = year)) +
geom_bar(color = "black", stat = "identity",
position = position_dodge2(preserve = "single"), width = 20) +
theme_bw() +
labs(title = '1998-2022',
x = "Mean CPUE/hr",
y = "Gulf stream position anomaly")
Correlations and GAM
#correlations gsi and cpue
lm_gsi<-lm(GSI ~ mean_cpue, data=gsi.sfob)
summary(lm_gsi)##
## Call:
## lm(formula = GSI ~ mean_cpue, data = gsi.sfob)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.13994 -0.50126 0.09577 0.57067 1.28856
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.80239 0.07705 10.413 <0.0000000000000002 ***
## mean_cpue -0.04394 0.01829 -2.403 0.0175 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7425 on 145 degrees of freedom
## Multiple R-squared: 0.03829, Adjusted R-squared: 0.03166
## F-statistic: 5.773 on 1 and 145 DF, p-value: 0.01754plot(GSI ~ mean_cpue, data=gsi.sfob, pch=1, col="dodgerblue") + abline(lm_gsi)
## integer(0)cor.test(gsi.sfob$mean_cpue, gsi.sfob$GSI, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: gsi.sfob$mean_cpue and gsi.sfob$GSI
## t = -2.4028, df = 145, p-value = 0.01754
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.34659356 -0.03489219
## sample estimates:
## cor
## -0.1956801ggscatter(gsi.sfob, x = "mean_cpue", y = "GSI",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.coeff.args = list(method="pearson", label.x=10, label.y=-1.0),
add.params=list(color="dodgerblue", fill="lightgray"),
xlab = "Mean CPUE/hr", ylab = "Mean Monthly GSI")
#GSI GAM
gsi_gam <- gam(mean_cpue ~ s(GSI), data=gsi.sfob, method="REML")
#gam.check(gsi_gam)
#summary(gsi_gam)
draw(gsi_gam, residuals=TRUE)
Food availablity
Larval tilefish eat zooplankton, likely calanus. Calanus finmarchicus are a copepod (crustacean) with a one-year life cycle and are an important food source for many commercially important species. Calanus spp. are lipid rich, herbivorous species that eat phytoplankton, diatoms in particular (Hobbs et al. 2020).
Diatoms are often represented as microplankton (>20 µm), but many species are of the nanoplankton size class (2-20 µm), and a smaller few may even overlap with picoplanton size class (<2 µm).
Calanus
Calanus is not as common in MAB, need to figure out dominant zooplankton in MAB.
# Calanus
calanus <- ecodata::calanus_stage %>% filter(Time %in% c(1998:2021))%>%
rename_all(., .funs = tolower) %>%
mutate(year = time)
ggplot() +
geom_line(data = calanus %>% filter(epu == 'GB',
var == 'adt-Spring'),
aes(x = year , y = value, col = epu), lwd = 1) +
geom_line(data = calanus %>% filter(epu == 'MAB',
var == 'adt-Spring'),
aes(x = year , y = value, col = epu), lwd = 1) +
labs(color = c('EPU')) +
theme_minimal()
Georges Bank
# GB
c5.gb <- calanus %>% filter(epu == 'GB', var == 'c5-Spring')
adult.gb <- calanus %>% filter(epu == 'GB', var == 'adt-Spring' )
df.c5 <- dplyr::full_join(recruit, c5.gb, by = join_by(year)) %>%
dplyr::select(year, recruit_est, value) %>%
tidyr::drop_na()
df.adt <- dplyr::full_join(recruit, adult.gb, by = join_by(year)) %>%
dplyr::select(year, recruit_est, value) %>%
tidyr::drop_na()
# Regression
ggscatter(df.c5, x = 'recruit_est', y = 'value',
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Recruitment estimate",
ylab = "Calanus c5 spring (No. per 100m^-3)",
title = 'c5') 
ggscatter(df.adt, x = 'recruit_est', y = 'value',
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Recruitment estimate",
ylab = "Calanus adult spring (No. per 100m^-3)",
title = 'Adult') 
# GLM
eqn <- as.formula(paste('recruit_est ~', paste(colnames(df.c5)[1],
collapse = " + ")))
mod0 <- glm(recruit_est ~ 1,
data = df.c5,
family = "poisson")
mod1 <- glm(eqn,
data = df.c5,
family = "poisson")
summary(mod0)##
## Call:
## glm(formula = recruit_est ~ 1, family = "poisson", data = df.c5)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 14.1979354 0.0001802 78775 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 8933474 on 20 degrees of freedom
## Residual deviance: 8933474 on 20 degrees of freedom
## AIC: 8933809
##
## Number of Fisher Scoring iterations: 4summary(mod1)##
## Call:
## glm(formula = eqn, family = "poisson", data = df.c5)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 43.50429247 0.05804867 749.4 <0.0000000000000002 ***
## year -0.01459584 0.00002891 -504.8 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 8933474 on 20 degrees of freedom
## Residual deviance: 8677368 on 19 degrees of freedom
## AIC: 8677705
##
## Number of Fisher Scoring iterations: 4AIC(mod0, mod1) %>% dplyr::arrange(AIC)## df AIC
## mod1 2 8677705
## mod0 1 8933809null_prediction <- exp(predict(mod0))
mod_prediction <- exp(predict(mod1))
plot(df.c5$year, df.c5$recruit_est, type = 'l')
lines(df.c5$year, null_prediction, col = "gray")
lines(df.c5$year, mod_prediction, col = "red")
Mid-atlantic
# Mid-Atlantic Bight
c5.mab <- calanus %>% filter(epu == 'MAB', var == 'c5-Spring')
adult.mab <- calanus %>% filter(epu == 'MAB', var == 'adt-Spring' )
df.c5 <- dplyr::full_join(recruit, c5.mab, by = join_by(year)) %>%
dplyr::select(year, recruit_est, value) %>%
tidyr::drop_na()
df.adt <- dplyr::full_join(recruit, adult.mab, by = join_by(year)) %>%
dplyr::select(year, recruit_est, value) %>%
tidyr::drop_na()
# Regression
ggscatter(df.c5, x = 'recruit_est', y = 'value',
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Recruitment estimate",
ylab = "Calanus c5 spring (No. per 100m^-3)",
title = 'c5') 
ggscatter(df.adt, x = 'recruit_est', y = 'value',
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Recruitment estimate",
ylab = "Calanus adult spring (No. per 100m^-3)",
title = 'Adult') 
Microplankton
# microplankton
micro<-read.csv(here::here('data/phyto_size_class/microplankton_ts_gtf_strata.csv'))Microplankton maps
Microplankton analysis
# microplankton by month with lag
micro.lag <- micro %>%
mutate(mean_micro_lag2 = lag(weighted_mean_micro,2),
mean_micro_lag3 = lag(weighted_mean_micro,3),
mean_micro_lag6 = lag(weighted_mean_micro,6))
# filter sfob.env
sfob <- sfob.env %>% filter(depth > 50, cpue_hr > 0) %>%
group_by(year,month) %>%
summarise(mean_cpue = mean(cpue_hr))
# Join with sf/ob data
micro.sfob <- dplyr::full_join(micro.lag, sfob, by = join_by(year, month)) %>%
dplyr::select(year, month, mean_cpue, mean_micro, weighted_mean_micro, mean_micro_lag2, mean_micro_lag3, mean_micro_lag6) %>%
tidyr::drop_na()
## Microplankton no lag
ggplot2::ggplot(micro.sfob,
aes(x=mean_cpue, y=weighted_mean_micro)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean Microplankton') +
labs(title = 'Mean Microplankton no lag') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Microplankton summer
ggplot2::ggplot(micro.sfob %>% filter(month %in% c(7,8,9)),
aes(x=mean_cpue, y=weighted_mean_micro)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean Microplankton') +
labs(title = 'Summer Mean Microplankton') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Microplankton spring
ggplot2::ggplot(micro.sfob %>% filter(month %in% c(4,5,6)),
aes(x=mean_cpue, y=weighted_mean_micro)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean Microplankton') +
labs(title = 'Spring Mean Microplankton') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## 2 month lag
ggplot2::ggplot(micro.sfob,
aes(x=mean_cpue, y=mean_micro_lag2)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean Microplankton') +
labs(title = 'Mean Microplankton 2 month lag') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## 3 month lag
ggplot2::ggplot(micro.sfob,
aes(x=mean_cpue, y=mean_micro_lag3)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean Microplankton') +
labs(title = 'Mean Microplankton 3 month lag') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## 6 month lag
ggplot2::ggplot(micro.sfob,
aes(x=mean_cpue, y=mean_micro_lag6)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean Microplankton') +
labs(title = 'Mean Microplankton 6 month lag') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
Micro correlation and GAM
# Correlation and plot
lm_micro<-lm(weighted_mean_micro ~ mean_cpue, data=micro.sfob)
summary(lm_micro)##
## Call:
## lm(formula = weighted_mean_micro ~ mean_cpue, data = micro.sfob)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.32848 -0.11240 -0.01818 0.07787 0.79427
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.270565 0.019885 13.607 < 0.0000000000000002 ***
## mean_cpue 0.007265 0.002372 3.062 0.00256 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1695 on 166 degrees of freedom
## Multiple R-squared: 0.05347, Adjusted R-squared: 0.04777
## F-statistic: 9.378 on 1 and 166 DF, p-value: 0.002563plot(weighted_mean_micro ~ mean_cpue, data=micro.sfob, pch=1, col="dodgerblue") + abline(lm_micro)
## integer(0)cor.test(micro.sfob$mean_cpue, micro.sfob$weighted_mean_micro, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: micro.sfob$mean_cpue and micro.sfob$weighted_mean_micro
## t = 3.0623, df = 166, p-value = 0.002563
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.08272836 0.36970715
## sample estimates:
## cor
## 0.2312416ggscatter(micro.sfob, x = "mean_cpue", y = "weighted_mean_micro",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.coeff.args = list(method="pearson"),
add.params=list(color="dodgerblue", fill="lightgray"),
xlab = "Mean CPUE/hr", ylab = "Weighted Mean Microplankton")
# GAM
micro_gam <- gam(mean_cpue ~ s(weighted_mean_micro), data=micro.sfob, method="REML")
#gam.check(micro_gam)
#summary(micro_gam)
draw(micro_gam, residuals=TRUE)
Chlorophyll-A
# CHL-A
chl<-read.csv(here::here('data/chl/chl_ts_gtf_strata.csv'))Maps (Chlorophyll-A)
Chlorophyll-A analysis
# chl by month with lag
chl.lag <- chl %>%
mutate(mean_chl_lag2 = lag(weighted_mean_chl,2),
mean_chl_lag3 = lag(weighted_mean_chl,3),
mean_chl_lag6 = lag(weighted_mean_chl,6))
# filter sfob.env
sfob <- sfob.env %>% filter(depth > 50, cpue_hr > 0) %>%
group_by(year,month) %>%
summarise(mean_cpue = mean(cpue_hr))
# Join with sf/ob data
chl.sfob <- dplyr::full_join(chl.lag, sfob, by = join_by(year, month)) %>%
dplyr::select(year, month, mean_cpue, mean_chl, weighted_mean_chl, mean_chl_lag2, mean_chl_lag3, mean_chl_lag6) %>%
tidyr::drop_na()
## CHL-a no lag
ggplot2::ggplot(chl.sfob,
aes(x=mean_cpue, y=weighted_mean_chl)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean CHL-a') +
labs(title = 'Mean CHL-a no lag') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## CHL summer
ggplot2::ggplot(chl.sfob %>% filter(month %in% c(7,8,9)),
aes(x=mean_cpue, y=weighted_mean_chl)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean CHL-a') +
labs(title = 'Summer Mean CHL-a') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## CHL fall
ggplot2::ggplot(chl.sfob %>% filter(month %in% c(10,11,12)),
aes(x=mean_cpue, y=weighted_mean_chl)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean CHL-a') +
labs(title = 'Fall Mean CHL-a') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## 2 month lag
ggplot2::ggplot(chl.sfob,
aes(x=mean_cpue, y=mean_chl_lag2)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean CHL-a') +
labs(title = 'Mean CHL-a 2 month lag') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## 3 month lag
ggplot2::ggplot(chl.sfob,
aes(x=mean_cpue, y=mean_chl_lag3)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean CHL-a') +
labs(title = 'Mean CHL-a 3 month lag') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## 6 month lag
ggplot2::ggplot(chl.sfob,
aes(x=mean_cpue, y=mean_chl_lag6)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean CHL-a') +
labs(title = 'Mean CHL-a 6 month lag') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
CHL-a and CPUE correlation
lm_chl<-lm(weighted_mean_chl ~ mean_cpue, data=chl.sfob)
summary(lm_chl)##
## Call:
## lm(formula = weighted_mean_chl ~ mean_cpue, data = chl.sfob)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.41118 -0.13268 -0.02328 0.11581 1.17495
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.625866 0.027749 22.554 <0.0000000000000002 ***
## mean_cpue 0.007907 0.003311 2.388 0.018 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2366 on 166 degrees of freedom
## Multiple R-squared: 0.03322, Adjusted R-squared: 0.0274
## F-statistic: 5.705 on 1 and 166 DF, p-value: 0.01804plot(weighted_mean_chl ~ mean_cpue, data=chl.sfob, pch=1, col="dodgerblue") + abline(lm_chl)
## integer(0)cor.test(chl.sfob$mean_cpue, chl.sfob$weighted_mean_chl, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: chl.sfob$mean_cpue and chl.sfob$weighted_mean_chl
## t = 2.3884, df = 166, p-value = 0.01804
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.0317402 0.3247221
## sample estimates:
## cor
## 0.182274ggscatter(chl.sfob, x = "mean_cpue", y = "weighted_mean_chl",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.coeff.args = list(method="pearson"),
add.params=list(color="dodgerblue", fill="lightgray"),
xlab = "Mean CPUE/hr", ylab = "Mean Monthly CHL-a")
CHL-a GAM
chl_gam <- gam(mean_cpue ~ s(weighted_mean_chl), data=chl.sfob, method="REML")
#gam.check(chl_gam)
#summary(chl_gam)
draw(chl_gam, residuals=TRUE)
SST Fronts
# SST fronts
fprob<-read.csv(here::here('data/sst_fronts/fprob_seasonal_ts_gtf_3x3.csv')) #fprob across all individual strata
fprob.ind<-read.csv(here::here('data/sst_fronts/fprob_seasonal_ts_gtf_indv_substrata_3x3.csv')) #mean for each of 14 substrata
fprob.ns<-read.csv(here::here('data/sst_fronts/fprob_seasonal_ts_gtf_n_v_s_substrata_3x3.csv')) #mean for substrata N and S of Hudson canyonMaps (SST fronts)
SST Fronts analysis
FPROB across all individual strata
#sfob by season
sfob_season<-read.csv(here::here('data/sst_fronts/sfob_byseason.csv'))
#join with fprob
fprob.sfob <- dplyr::full_join(fprob, sfob_season, by = join_by(year, season)) %>%
dplyr::select(year, season, mean_cpue, mean_fprob, weighted_mean_fprob) %>%
tidyr::drop_na()
## All months
ggplot2::ggplot(fprob.sfob,
aes(x=mean_cpue, y=weighted_mean_fprob)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean SST frontal probability') +
labs(title = 'Mean SST frontal probability') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Winter
ggplot2::ggplot(fprob.sfob %>% filter(season == "winter") ,
aes(x=mean_cpue, y=weighted_mean_fprob)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean SST frontal probability') +
labs(title = 'Winter Mean SST frontal probability') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Spring
ggplot2::ggplot(fprob.sfob %>% filter(season == "spring"),
aes(x=mean_cpue, y=weighted_mean_fprob)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean SST frontal probability') +
labs(title = 'Spring Mean SST frontal probability') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Summer
ggplot2::ggplot(fprob.sfob %>% filter(season == "summer"),
aes(x=mean_cpue, y=weighted_mean_fprob)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean SST frontal probability') +
labs(title = 'Summer Mean SST frontal probability') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Fall
ggplot2::ggplot(fprob.sfob %>% filter(season == "fall"),
aes(x=mean_cpue, y=weighted_mean_fprob)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean SST frontal probability') +
labs(title = 'Fall Mean SST frontal probability') +
stat_cor(aes(label=..rr.label..)) +
theme_bw()
## Correlation and GAM
lm_fprob<-lm(mean_cpue ~ weighted_mean_fprob, data=fprob.sfob)
summary(lm_fprob)##
## Call:
## lm(formula = mean_cpue ~ weighted_mean_fprob, data = fprob.sfob)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.089 -3.557 -1.112 1.547 26.933
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 18.372 9.781 1.878 0.0642 .
## weighted_mean_fprob -14.933 12.527 -1.192 0.2369
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.586 on 76 degrees of freedom
## Multiple R-squared: 0.01835, Adjusted R-squared: 0.005438
## F-statistic: 1.421 on 1 and 76 DF, p-value: 0.2369plot(weighted_mean_fprob ~ mean_cpue, data=fprob.sfob, pch=1, col="dodgerblue") + abline(lm_fprob)
## integer(0)cor.test(fprob.sfob$mean_cpue, fprob.sfob$weighted_mean_fprob, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: fprob.sfob$mean_cpue and fprob.sfob$weighted_mean_fprob
## t = -1.1921, df = 76, p-value = 0.2369
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.34753296 0.08975701
## sample estimates:
## cor
## -0.1354798ggscatter(fprob.sfob, x = "mean_cpue", y = "weighted_mean_fprob",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.coeff.args = list(method="pearson"),
add.params=list(color="dodgerblue", fill="lightgray"),
xlab = "Mean CPUE/hr", ylab = "Mean Seasonal SST Frontal probability")
fprob_gam <- gam(mean_cpue ~ s(weighted_mean_fprob), data=fprob.sfob, method="REML")
#gam.check(fprob_gam)
#summary(fprob_gam)
draw(fprob_gam, residuals=TRUE)
FPROB of each substrata
#join with sfob
ind.sfob <- dplyr::full_join(fprob.ind, sfob_season, by = join_by(year, season)) %>%
dplyr::select(year, season, id, mean_cpue, mean_fprob, weighted_mean_fprob) %>%
tidyr::drop_na()
# fprob by substrata
ggplot2::ggplot(ind.sfob,
aes(x=mean_cpue, y=weighted_mean_fprob)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean SST frontal probability') +
labs(title = 'Mean SST frontal probability') +
stat_poly_eq(aes(label = paste(..rr.label.., sep = "~~~")),
label.x.npc = 30, label.y.npc = 0.15,
parse = TRUE, size = 3) +
theme_bw() +
facet_wrap(~id) +
theme_facet()
## Winter
ggplot2::ggplot(ind.sfob %>% filter(season == "winter") ,
aes(x=mean_cpue, y=weighted_mean_fprob)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean SST frontal probability') +
labs(title = 'Winter Mean SST frontal probability') +
stat_poly_eq(aes(label = paste(..rr.label.., sep = "~~~")),
label.x.npc = 30, label.y.npc = 0.15,
parse = TRUE, size = 3) +
theme_bw() +
facet_wrap(~id) +
theme_facet()
## Spring
ggplot2::ggplot(ind.sfob %>% filter(season == "spring") ,
aes(x=mean_cpue, y=weighted_mean_fprob)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean SST frontal probability') +
labs(title = 'Spring Mean SST frontal probability') +
stat_poly_eq(aes(label = paste(..rr.label.., sep = "~~~")),
label.x.npc = 30, label.y.npc = 0.15,
parse = TRUE, size = 3) +
theme_bw() +
facet_wrap(~id) +
theme_facet()
## Summer
ggplot2::ggplot(ind.sfob %>% filter(season == "summer") ,
aes(x=mean_cpue, y=weighted_mean_fprob)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean SST frontal probability') +
labs(title = 'Summer Mean SST frontal probability') +
stat_poly_eq(aes(label = paste(..rr.label.., sep = "~~~")),
label.x.npc = 30, label.y.npc = 0.15,
parse = TRUE, size = 3) +
theme_bw() +
facet_wrap(~id) +
theme_facet()
## Fall
ggplot2::ggplot(ind.sfob %>% filter(season == "fall") ,
aes(x=mean_cpue, y=weighted_mean_fprob)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean SST frontal probability') +
labs(title = 'Fall Mean SST frontal probability') +
stat_poly_eq(aes(label = paste(..rr.label.., sep = "~~~")),
label.x.npc = 30, label.y.npc = 0.15,
parse = TRUE, size = 3) +
theme_bw() +
facet_wrap(~id) +
theme_facet()
## Scatter plot with correlations
#plot cuts off p values unless you open it in a new window and save
ggscatter(ind.sfob, x = "mean_cpue", y = "weighted_mean_fprob",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.coeff.args = list(method="pearson",
label.x=20, label.y=0.6, size =3),
add.params=list(color="dodgerblue", fill="lightgray"),
xlab = "Mean CPUE/hr", ylab = "Mean SST frontal probability") +
facet_wrap(~id) +
theme_facet()
N/S substrata
fprob_ns<-read.csv(here::here('data/sst_fronts/fprob_ns.csv'))
fprob.ns<-fprob_ns[!(is.na(fprob_ns$mean_cpue)), ]
# fprob by N/S Hudson Canyon
ggplot2::ggplot(fprob.ns,
aes(x=mean_cpue, y=weighted_mean_fprob)) +
geom_point(color = 'black') +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
xlab('Mean CPUE/hr') +
ylab('Weighted mean SST frontal probability') +
labs(title = 'Mean SST frontal probability') +
stat_poly_eq(aes(label = paste(..rr.label.., sep = "~~~")),
label.x.npc = 10,
parse = TRUE, size = 5) +
theme_bw() +
facet_wrap(~substrat) +
theme_facet()
## Scatter plot with correlations
ggscatter(fprob.ns, x = "mean_cpue", y = "weighted_mean_fprob",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.coeff.args = list(method="pearson",
label.x=20, label.y=0.6, size =3),
add.params=list(color="dodgerblue", fill="lightgray"),
xlab = "Mean CPUE/hr", ylab = "Mean SST frontal probability") +
facet_wrap(~substrat) +
theme_facet()
References:
Hobbs, L., Banas, N. S., Cottier, F. R., Berge, J., & Daase, M. (2020). Eat or sleep: availability of winter prey explains mid-winter and spring activity in an Arctic Calanus population. Frontiers in Marine Science, 7, 541564.